// Mesh size
//h = 8.0; // implies 1846 elements
//h = 4.0; // implies 7826 elements (optimize Netgen)
//h = 2.0; // implies 36022 elements

h = 1.0;

h_cil = h/8;

// Define the cylinder
xcn = 0.0;
ycn = 0.0;
rad = 0.5;
zbs = 3.14*(2*rad);
// the cylinder is defined using four arcs
cen = newp; Point(cen) = {xcn,ycn,zbs,h};
A = newp; Point(A) = {xcn,ycn-rad,zbs,h_cil};
B = newp; Point(B) = {xcn+rad,ycn,zbs,h_cil};
C = newp; Point(C) = {xcn,ycn+rad,zbs,h_cil};
D = newp; Point(D) = {xcn-rad,ycn,zbs,h_cil};

a1 = newreg; Circle(a1) = {A,cen,B};
a2 = newreg; Circle(a2) = {B,cen,C};
a3 = newreg; Circle(a3) = {C,cen,D};
a4 = newreg; Circle(a4) = {D,cen,A};

hol = newreg; Line Loop(hol) = {-a4,-a3,-a2,-a1}; 

// Now the outer box
dxm =  5.0;
dxp = 20.0;
dy  =  7.0;
cbl = newp; Point(cbl) = {-dxm,-dy,zbs,h};
cbr = newp; Point(cbr) = { dxp,-dy,zbs,h};
ctr = newp; Point(ctr) = { dxp, dy,zbs,h};
ctl = newp; Point(ctl) = {-dxm, dy,zbs,h};

bot = newc; Line(bot) = {cbl,cbr}; 
rig = newc; Line(rig) = {cbr,ctr}; 
top = newc; Line(top) = {ctr,ctl}; 
lef = newc; Line(lef) = {ctl,cbl}; 

box = newreg; Line Loop(box) = {bot,rig,top,lef}; 

// Define the 2D domain
domain2D = news; Plane Surface(domain2D) = {box,hol};

// Extrusion to obtain the 3D domain
// spanwise direction 2*pi*D as in Snyder-Degrez-2003
domain3D = news; Extrude{0,0,2*3.14*(2*rad)}{Surface{domain2D};}

// Central refinement
//Field[1] = MathEval;
//Field[1].F = Sprintf( "min( %g/2*(0.25 + 1.75*(y/%g)^2 + (x/%g)^4) , %g )" , h , dy , dxp , (2*3.14*(2*rad))/2.5 );
//Field[1].F = Sprintf( "0.3*(1+abs(x)/5)*(abs(y)/3+0.5)");
//Field[1].F = Sprintf( "1.5*(0.1*abs(z)+0.0001*4*abs(z)/3.14+0.4)");
//Field[1].F = Sprintf( "(abs(z)/3.14+0.4)*%g/2*(0.25 + 1.75*(y/%g)^2 + (x/%g)^4)" , h , dy , dxp );
//Background Field = 1;

// Merge a post-processing view containing the target mesh sizes
Merge "abc.pos";

// Apply the view as the current background mesh
Background Mesh View[0];
